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Preface 


Since  the  early  1960 's,  science  has  made  great  strides  in  the 
development  of  high  powered  laser  systems  and  in  the  theory  of  optimal 
estimation  and  control.  This  report  attempts  to  combine  these  areas 
of  technology  by  presenting  an  extended  Kalman  filter  that  cem  be 
used,  together  with  the  outputs  of  a forward  looking  infrared  sensor, 
in  an  open-loop  tracking  problem.  This  subject  was  chosen  because  of 
the  great  significance  of  the  high  energy  laser  as  a potenticJ.  weapon 
system  and  because  of  the  need  for  research  in  this  highly  technical 
area.  Hopefully,  the  topic  presented  here  will  benefit  the  personnel 
at  the  Air  Force  weapons  Laboratory  and  aid  them  in  designing  a better 
weapons  system. 

My  thanks  go  out  to  all  those  people  whose  cissistance  and  guidcuice 
was  given  to  me  and  without  whose  help  this  would  have  been  in\possible. 
Mr.  Peter  Maybeck  was  exceptionally  helpful  as  I(^y  advisor.  His  meUiy 
comments  were  beneficial  and  his  door  was  always  open  for  advice. 
Additional  gratitude  must  be  expressed  to  Captains  Stan  Robinson,  Gary 
Reid,  and  Willieun  Wiesel  for  their  interest  and  support  in  conipleting 
this  thesis . A sincere  thax>ks  also  goes  to  Captains  stem  Lewantowlcz 
and  Felix  Morgan  at  the  Air  Force  Weapons  Laboratory  for  their  helpful 
suggestions  euid  for  the  time  they  spent  digging  up  letters  2md  reports 
pertenent  to  this  study. 

The  final  debt  of  thcUiks  I owe  is  to  my  wife  Diane  for  her  in- 
credible undcrsteuiding,  her  enormous  love,  her  gentle  prodding,  and 
for  the  many  hours  she  spent  alone  during  this  period.  v;e  were  in  this 
together . 

Daniel  e.  Mercier 
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Abstract 

i 

An  extended  Kalman  filter  algorithm,  using  outputs  from  a for- 
ward looJd-ng  infrared  (rEI?'  sensor  as  measurements,  is  used  to  track 
a point  source  target  in  an  open  loop  tracking  problem.  The  filter 
estimates  the  translational  position  changes  of  the  target  in  the  FLIR 
field  of  view  due  to  two  effects;  actual  target  motion,  and  apparent 
motion  caused  by  atmospheric  turbulence . Sixteen  cases  are  examined 
to  determine  the  performance  of  the  filter  as  a function  of  signal-to- 
noise  ratio,  Gaussian  beam  size,  the  ratio  of  rms  target  motion  to 
RMS  atmospheric  jitter,  target  correlation  times,  and  mismatches 
between  the  true  shape  and  the  shape  assumed  by  the  filter.  The  per- 
formance of  the  extended  Kalman  filter  is  compared  to  the  performance 
of  an  existing  correlation  tracker  under  identical  initial  conditions. 

^ A one  sigma  tracking  error  of  .2  and  .8  picture  elements  is  obtained 

with  signal-to-noise  ratios  of  20  and  1 respectively.  No  degradation 
in  performanje  is  observed  when  the  beam  size  is  decreased  or  when  the 
target  correlation  time  is  increased  over  a limited  range,  when  filter 
parameters  are  adjusted  to  reflect  this  knowledge.  Sensitivity  analysis 
shows  that  the  filter  is  robust  to  minor  changes  in  target  intensity 
size. 
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an  extended  KALMAN  FILTER  FOR 

I 

USE  IN  A SHARED  APERTURE 
MEDIUM  RANGE  TRACKER 

I.  Introduction 

Background 

Since  the  early  1960 's,  there  has  been  an  impressive  growth  in 
laser  technology  throughout  the  world.  Now,  only  fifteen  years  later, 
the  laser  beam  has  become  a prime  candidate  as  a potential  weapon  system 
and  is  currently  being  tested  against  airborne  targets . As  with  all 
new  systems,  however,  many  problems  must  be  resolved  before  the  high 
^ energy  laser  can  be  used  as  em  effective  weapon. 

Problem 

One  of  these  problems  is  the  precision  pointing  of  the  laser  at  a 
given  target.  The  Air  Force  Weapons  Laboratory  is  currently  examining 
several  methods  of  tracking  a target  for  the  purpose  of  depositing  high 
power  laser  energy  on  the  target  in  the  presence  of  several  disturbances . 

These  disturbances  include  any  effect  that  can  cause  relative  motion  ' 

between  the  beam  and  the  target,  such  as  target  motion,  mirror  vibration, 
and  atmospheric  jitter. 

One  tracking  method  currently  being  used  employs  a forward  looking 
infrared  sensor  (FLIR)  together  with  a correlation  algorithm  to  provide 
relative  target  position  information  to  the  laser  pointing  system.  The 
l'  algorithm  provides  this  information  by  first  storing  a complete  set  of 

1 


target  intensity  data  from  the  outputs  of  the  FLIR,  and  then  correlating 
that  data  with  new  information  at  a later  time.  Through  the  use  of 
cross  correlation  techniques,  the  algorithm  can  estimate  the  relative 
position  offset  from  one  set  of  data  to  the  next.  Using  this  relative 
position  information,  commands  can  then  be  sent  to  the  gimbals  controlling 
the  line  of  sight  of  the  infrared  system  to  keep  it  centered  on  the  target. 
This  centering  of  the  FLIR  on  the  target  also  physically  points  the  laser 
towards  the  target  since  the  high  energy  laser  is  servoed  to  the  FLIR 
system.  This  type  of  tracker  needs  no  prior  information  about  the  type 
of  target  to  perform  the  tracking  function,  and  therefore,  it  is  well 
suited  to  many  general  applications. 

In  many  practical  tracking  problems,  however,  the  type  of  target 
being  tracked  will  be  known,  if  even  in  a very  general  sense.  This  im- 
plies that  certain  target  parameters  such  as  shape , size , and  accelera- 
tion characteristics,  will  either  be  known  or  could  be  estimated.  Also, 
the  statistical  effects  of  atmospheric  distortion  on  radiated  wavefronts 
are  known  and  could  supply  information  to  a tracker  that  would  aid  in 
separating  the  true  target  relative  motion  from  the  apparent  motion  due 
to  the  atmospheric  disturbance . This  separation  is  important  since  the 
wavefront  of  the  high  energy  laser  will  not  undergo  the  same  distortion 
in  the  atmosphere  as  the  infrared  wavefronts  emamating  from  the  target. 

This  readily  available  information,  not  used  by  the  correlation  tracker, 
could  be  used,  together  with  the  outputs  of  the  FLIR  system,  in  an 
algorithm  which  would  optimally  estimate  relative  target  position  infor- 
mation. One  way  to  do  this  would  be  to  use  an  extended  Kalman  filter 
in  the  tracking  loop  in  place  of  the  present  correlation  algorithm. 
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Research  Goals 


This  report  presents  the  use  of  an  extended  Kcilman  filter  in  an 
open  loop  tracking  problem  for  the  eventual  use  in  the  precision  pointing 
of  a laser  beam  at  a target.  The  study  was  performed  by  simulating  the 
filter  algorithm  and  the  PLIR  outputs  on  a CDC-6600  computer.  A variety 
of  simulations  were  performed,  examining  scenarios  with  varying  signal- 
to-noise  ratios  (S/N) , target  intensity  pattern  sizes,  and  ratios  of 
target  dynamics  to  atmospheric  jitter. 

As  a basis  for  comparison,  a correlation  tracking  algorithm,  de- 
vised by  the  Aeronutronic  Ford  Corporation,  using  single-degree-of-free- 
dom  correlation  techniques  in  two  independent  directions , was  also  sim- 
ulated. Both  algorithms  were  r\in  with  identical  inputs  eind  with  iden- 
tical information  from  the  FLIR.  The  extended  Kalman  filtering  algorithm, 
however,  made  use  of  some  apriori  knowledge  of  the  size,  shape,  and 
typical  dynamics  of  the  target  along  with  some  knowledge  of  the  jitter 
characteristics  of  the  atmosphere . Outputs  for  the  two  algorithms  were 
then  compared  to  evaluate  the  performance  of  each  under  varying  condi- 
tions . Several  simulations  were  also  made  to  investigate  the  sensitivity 
of  the  extended  Kalman  filter  algorithm  to  incorrectly  made  assumptions . 

In  order  to  limit  the  problem  to  one  which  could  be  accoirplished  in  a 
specified  allotment  of  time,  only  the  open  loop  tracking  problem  was 
addressed.  In  this  way,  controller  design  could  be  ler  as  a follow-on 
study . 

Assumptions 

Target . Most  applications  for  which  the  laser  system  is  being  con- 
sidered require  the  acquisition  and  tracking  of  targets  at  long  r2Uiges. 
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Because  of  these  distances,  even  very  large  targets  appear  as  point 
sources  of  infrared  radiation  and  can  be  accurately  modelled  as  such. 

This  was  the  only  type  of  target  considered  since  a variation  in  target 
shape  would  require  a major  rewrite  of  those  subroutines  simulating 
the  FLIR  outputs . 

FLIR  Sensor . The  outputs  of  a typical  forward  looking  infrared 
sensor  are  continuous  time  averages  cf  an  array  of  infrared  detectors 
as  they  are  mechanically  scanned  through  a limited  field  of  view  (FOV) . 
As  the  detectors  are  scanned  in  both  the  horizontal  and  vertical  direc- 
tions, they  emit  electrical  current  proportional  to  the  number  of  photons 
entering  the  face  of  the  detectors.  In  order  to  get  an  accurate  repre- 
sentation of  the  average  number  of  photons  present  at  a specific  point 
in  space , the  currents  from  each  detector  as  it  passes  that  point  are 
combined,  electronically  averaged,  and  converted  to  digital  signals  for 
processing.  This  electronic  averaging  is  performed  in  real  time  by 
processing  the  output  of  each  detector  through  a delay  line  and  into  a 
summing  junction.  This  sximmation  is  then  divided  by  the  number  of  de- 
tectors in  the  array  to  obtain  the  final  averaged  output.  This  averaging 
technique,  as  shown  in  Figure  1,  compensates  for  the  time  difference 
that  each  detector  is  scemned  across  a specific  point  in  space.  The  out- 
put of  the  averaging  process  can  then  either  be  stored  or  displayed  on 
a cathode  ray  tube  (CRT)  in  real  time.  Each  output  will  correspond  to 
one  picture  element  on  the  CRT,  usually  termed  pixel.  The  horizontal 
and  vertical  scanning  of  the  detectors  through  the  FLIR  field  of  view 
will  result  in  an  array  of  pixels  called  a frame  of  data.  Normal  frame 
rates  are  on  the  order  of  30  Hz.  For  the  purposes  of  this  study,  one 
frame  of  (fbta  will  be  considered  as  one  measurement  array,  llie  scan 
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Figure  1.  Electronic  Averaging  of  Detector  Outputs 


time,  therefore,  predetermines  the  upper  bound  on  the  measurement  rate. 
The  measurement  array  can  be  as  large  as  the  FLIR  field  of  view,  but 
because  of  the  fast  measurement  rates,  the  target  will  generally  not 
move  more  than  a few  pixels  between  measurements . This  permits  looking 
at  a small  pixel  array  for  tracking  purposes . In  this  study , an  8 by  8 
I array  of  picture  elements  was  used  as  the  tracking  window,  resulting  in 

a tolerable  amoxint  of  convputer  storage  area. 

Glint.  As  mentioned  above,  the  target  is  assumed  to  be  a point 
source  with  a time  invariemt  radiation  intensity  pattern,  normeJ-ly  termed 
glint.  Due  to  the  physics  of  wave  propagation,  this  glint  has  a Gaussi2ui 
distribution.  This  distribution  is  the  function  of  two  Euler  eUigles 
which  describe  the  orientation  of  the  boresight  of  the  FLIR  with  respect 
to  the  path  of  maximum  glint  emanating  from  tl»e  target.  The  Euler 
2mgles  are  the  actual  pointing  errors  in  the  horizontal  and  vertical 
directions.  These  pointing  errors  translate,  in  general,  to  two  linear 
pointing  errors  in  the  FLIR  field  of  view  as  shown  in  Figure  2. 

Hie  glint  reflected  from  the  target  can  now  be  modelled  as  a bi- 
^ variate  Gaussian  function  with  a two-dimensional  pointing  error  argu- 

ment. The  intensity  anywhere  in  the  nieasurement  plane  can  be  described 


I 

I 


as 


(1) 


n- 


l<x,y)  = [(x-Xp^^)^  + {y-ypeak)^]>  + *=(x,y) 


= the  maximum  target  glint 

= the  dispersion  of  the  Gaussieu:i  function 

= the  coordinates  cf  any  point  in  the  measurement 
array 

^^^ak'Vpeak^  “ coordinates  of  the  center  of  the  Gaussian 
intensity  function 

b(x,y)  = background  noise  as  a function  of  position  in 

the  measurement  plane 

Figure  3 shows  this  two-dimensional  glint  model  without  the  background 
noise  term.  Points  of  constant  intensity  for  this  ideal  function  appear 
as  circles  in  the  FLIR  field  of  view  as  shown  in  Figure  2. 

Algorithm.  Because  this  is  an  open  loop  tracking  problem,  the  real 
task  becomes  creating  an  algorithm  which  will  accurately  estimate  the 
position  of  the  point  of  maximum  intensity  with  respect  to  the  center 
of  the  8 by  8 array  of  pixels.  The  intensity  function  is,  of  course, 
moving  around  in  the  field  of  view  due  to  the  dynamics  of  the  teurget  emd 
also  due  to  the  atmospheric  jitter  which  is  introduced  as  the  radiated 
waves  propagate  through  the  atmosphere.  The  center  of  the  intensity 
function  or  the  FLIR  image  plane  is,  therefore,  described  by; 


where 


^max 

{x,y) 


^peak^^^ 


XpCt)  + x^(t) 
yo^t)  + y;^(t) 


(2) 
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where  Xp  and  are  due  to  target  dyneimics  and  X;^  and  y^^  are 

due  to  atmospheric  jitter.  An  ideal  tracking  algorithm  woula  accurately 
estimate  xp  and  yp  in  order  to  eliminate  the  effects  of  the  jitter. 

Vibrations  of  the  FLIR  system  could  also  cause  relative  pointing 
errors,  but  in  this  study,  the  FLIR  system  was  assumed  to  be  on  the  sur- 
face of  the  earth  or  some  other  stable  platform  so  that  the  vibration 
effects  were  neglected. 

Overview 

The  following  three  chapters  will  describe,  in  detail,  the  mathe- 
matical models  used  in  the  computer  simulation.  Chapters  II,  III,  and 
IV  will  present  the  truth  model,  filter  model,  and  correlation  algorithm 
respectively.  Chapter  V will  describe  the  computer  simulation  and  the 
error  analysis  techniques . Chapter  VI  will  show  the  results  cf  all 
the  test  cases.  Chapter  VII  will  present  the  conclusions  and  recommen- 
dations . 
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II.  Truth  McxJel 


Introduction 

The  problem,  as  stated  in  Chapter  1,  concerns  the  tracking  of  a 
moving  target  from  a stabilized  platform,  using  the  outputs  of  a FLIR 
sensor  as  measurements  of  target  radiation.  The  imaging  tracking  con- 
trol loop  for  the  high  energy  laser  has  a bandwidth  on  the  order  cf  2 
Hz.  This  control  loop  will,  therefore,  compensate  for  tracking  error 
only  in  this  range  of  interest.  Higher  order  aberrations  are  controlled 
separately  and  will  not  be  considered  here.  The  two  processes  causing 
movement  of  the  intensity  fvinction  in  the  low  frequency  range  are  target 
dynamics  and  atmospheric  jitter,  caused  by  the  constantly  changing  atmo- 
spheric conditions . 

Target  Dynamics 

The  purpose  of  any  tracking  algorithm  is  to  predict,  as  accurately 
as  possible,  the  movement  of  the  target  with  respect  to  the  center  of 
the  field  of  view  of  the  sensor.  Ideally,  when  the  control  loop  is 
closed,  the  target  should  remain  centered.  Imperfections  and  noise 
terras,  however,  will  combine  to  produce  errors  to  this  perfect  tracker. 

In  many  cases,  these  errors  can  be  well-described  as  outputs  of  first 
order  shaping  filters  driven  by  zero  mean,  white,  Gaussian  noise.  This 
type  of  model  tcikes  into  account  the  time-correlated  properties  of  the 
error  and  has  been  used  successfully  in  the  past  to  describe  this  type 
of  error.  Although  there  may  be  better  models  to  describe  the  errors 
of  specific  targets,  this  model  is  generally  applicable  and  simple  to 
implement  while  still  yielding  some  of  the  basic  characteristics  of 
trajectories . The  outputs  of  the  first  order  filters  can  be  expressed 
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mathematically  as 


3^(t)  = -l/T^Xj^(t)  + Xj^(t)  (3) 

yjj(t)  = -1/T^jj(t)  +W2(t)  (4) 

where 

E{wj^(t)}  «=  E{w2(t)}  = 0 (5) 

E{w2^(t)Wj^(s) } = E{w2(t)w2(s) } = 2a|/T^6(t-s)  (6) 


EiWj^(t)W2{s)  } = 0 

and 

Tiji  = correlation  time 

W2^(t),W2(t)  = independent,  white,  Gaussian  noise  processes 

= the  desired  variance  on  the  outputs  x^  and  y^ 

Atmospheric  Jitter  Model 

The  process  used  to  model  the  translational  position  changes  of 
the  intensity  distribution  due  to  disturbances  in  the  atmosphere,  was 
based  on  a study  by  the  Analytic  Sciences  Corporation  (TASC)  (Ref  1:29, 
30).  As  part  of  this  study,  TASC  generated  a power  spectral  density 
representation  for  the  effects  of  atmospheric  turbulence  using  a pro- 
gram supplied  by  the  Air  Force  Weapons  Laboratory  (Ref  2) , This  func- 
tion was  then  approximated  by  a third  order  shaping  filter  driven  by 
white,  Gaussiem  noise  to  arrive  at  a simple,  yet  adequate,  model. 
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Figure  4 presents  both  the  power  spectral  density  and  the  approxima- 
tion as  presented  in  the  TASC  report.  Using  this  same  approximation, 
the  target  intensity  jitter  in  each  independent  direction,  is  modeled 
as  the  output  of  the  filter  shown  in  Figure  5 . The  white  noise  process 
driving  the  filter,  was  given  the  following  characteristics; 


E{w(t)}  = 0 


(7) 


E{w(t)w(s)}  = 1 6(t-s) 


(8) 


The  value  of  gain,  K,  can  now  be  adjusted  to  obtain  a desired  root  mean 
squared  (RMS)  jitter  characteristic  on  the  output.  This  allows  studying 
the  effects  of  differing  amounts  of  jitter  corresponding  to  a wide 
range  of  atmospheric  conditions . A detailed  description  of  how  the 
actual  value  of  gain,  K,  is  determined  is  shown  in  Appendix  C. 


State  Space  Model 

Transferring  these  models  into  state  space  notation  yields  a time- 
invariant  vector  differential  equation  of  the  form 


x^(t)  = F^  x^(t)  + Gjj.  w^{t)  (9) 

where 

F “ the  truth  model  plant  matrix 
-T 

3^(t)  = the  truth  mot.  ,!,  state  vector 
= the  truth  model  input  matrix 
w„(t)  = a vector  of  white  Gaussian  noise  inputs 

and 


E{v^(t)}  - 0 


(10) 
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Figure  4.  Atmospheric  Turbulence  Spectral  Representation 


T 1 

E[v^(t)^(s)]  = 2.T'5(t-s) 


At  first  attempt,  a standard  phase  variable  form  of  the  equation  was 
used.  This  form,  however,  was  found  to  be  numerically  imprecise  in  the 
calculation  of  certain  desired  quantities  and  the  Jordon  canonical  form 
was  used  instead.  Detailed  information  concerning  the  numerical  problems 
encountered  can  be  found  in  Appendix  A. 

Using  the  Jordan  canonical  form,  as  developed  in  Appendix  B for  the 
atmospheric  turbulence  model,  equation  (9)  becomes 


2^(t) 


— 0000000 
T 

0 -a  000000 


00  -b  10000 


000  -b  0000 


0000  — 000 
^T 

00000  -a  00 


10  0 0 


0 0 0 

0 G2  0 0 


0 G^  0 0 


0 0 10 


0 0 0 G, 


Wj(t) 


000000  -bl 


0 0 0 G. 


0000000  -b 


0 0 0 G3  (13) 


where 


= KabV(a-b) 


“ "^1 


= Kab^/(a-b) 


and 


= x(l) 


= x(2)  + x(3) 


Yd  = x(5) 


= x(6)  + x(7) 


The  output  matrix  is  then 


y^Ct)  = 


Vak^t) 


ypeak^"^) 


11100000 

00001110 


(14) 


It  is  apparent  in  this  assumed  formulation,  that  the  horizontal  and 
vertical  channels  are  independent  with  identical  process  characteristics , 


Propagation  Equations 

The  solution  to  equation  (9)  is 


JCj(t)  = £3.(t,t^)}^(tjL)  + /^(t,T)Gj(T)^(T)  dT 

^i 


(15) 


where  ^^,(t,t^)  is  the  state  transition  matrix  and  satisfies  the  matrix 
differential  equation 


i(t,tj^)  = Fi^(t,ti) 


(16) 


and  where  ^(t^ft^)  = ^ (the  identity  matrix). 

Since  the  truth  model  plant  matrix  is  a constant,  tlie  state  transition 
matrix  is  a function  of  (t-tj^)  , and  therefore,  is  also  constant  for  a 
fixed  sampling  time . Because  of  the  independent  channels , the  state 
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transition  matrix  is  block  diagonal  with  both  4 by  4 blocks  equal  to 


-bAt  ^-bAt 


where  At  = (t^^+i  - t^)  , 

The  integral  term  shown  in  equation  (15)  is  a stochastic  integral, 
This  requires  the  probabilistic  characterization  of  the  integral  in 
order  to  describe  the  characteristics  of  the  process  >^(t)  . Knowing 
that  the  integral  has  a zero  mean  Gaussian  dis  -ribution,  and  letting 


W^  (ti)  = ^(T)  dT 

d tj_ 


it  can  be  determined  that  (Ref  4:174,181) 


G (T)  ^('r)G^{T)$'^(t.  ,T)dT 
— fT  i+X 


'V‘j> 


An  equivalent  discrete  time  model  to  portray  the  evolution  of  the 
2^(t)  process  can  new  be  given  as 


“ ^T2^T(ti)  + W (tj^) 


where 


k.-: . 


E[^^{ti)]  = 0 

:[^(ti)w^(tj)]  = 6.  . 


(23) 

(24) 


(note  that  = Qjp  (t^)  for  all  t^  because  of  stationary  inputs  and 

a time  invariant  system)  eind 


1+1  ip  ip 

^i 


(25) 


Another  way  of  writing  this  equation  is 


(26) 


where  ^ lower  triangular  matrix  known  as  the  Cholesky  square 

root  of  and  satisfies  the  equation 


■VSp.  " ■ V 


(27) 


and 


E|vj„(t.)n„(t.))  = 

Equation  (26)  has  the  same  properties  as  equation  (22)  since 

“ - 1^,1  v/2t^ 


(28) 


I’d  * -^d  - ■ “M 

Equation  (26)  propagates  the  truth  model  states  in  time  and  accounts 
for  the  stochastic  nature  of  the  problem. 


Measurements 


The  tjnith  model  measurements,  by  definition,  are  the  best 
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representations  possible  for  the  outputs  of  the  FLIR  sensor.  These  out- 


puts, as  discussed  earlier,  are  a measure  of  the  average  intensity  of 
tcirget  radiation  over  the  detector  surfaces  as  they  are  swept  across  the 
image  plane.  Sweep  rates  are  on  the  order  to  30  Hz  which  means  that 
measurements  are  available  at  discrete  time  points  every  1/30  seconds . 

The  input  radiation  consists  of  apparent  target  radiation  and 
FLIR  generated  noise  term  to  account  for  noisy  electronics , thermal 
radiation,  and  dark  current  noises  internal  to  the  detector  array. 

Given  all  this,  the  output  of  the  i-j-th  pixel  at  any  time  t , can  be 
described  as 


I 


Zij  (t) 


pixel 


+ n.  j (t)  + f . j (t) 


(30) 


where 


Zij  (t) 


(x,y) 


output  of  i-j-th  pixel  at  time  t 
area  of  pixel  under  consideration 
coordinates  of  any  point  in  the  pixel 


(^peak  ' 
ypeak^^")) 


coordinates  of  the  center  of  the  intensity 
distribution  with  respect  to  the  center  of  the 
FLIR  field  of  view  at  time  t (see  Eq.  (2)) 


nuoc 


”ij 


maximum  target  intensity 

dispersion  of  the  Gaussian  intensity  function 
background  noise  term  for  the  i-j-th  pixel 
FLIR  noise  term  for  the  i-j-th  pixel 
The  two  noise  terms,  (t)  and  f^ ^ (t)  , are,  in  general,  independent 

wide  band  noises  that  c^m  be  accurately  modeled,  according  to  the  Air 
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Force  Weapons  Laboratory,  as  zero  mean  white,  Gaussian  noise  terms  with 
the  following  statistics ; 


E[nkl(ti)]  = E 

|fki(ti)]  - 0 

for 

all  i,k,l 

(31) 

E[nki(ti)nki(tj)] 

= 

for 

all  i,k,l 

(32) 

E(fkl(ti)fki{tj) 

= a|6. . 

for 

all  i,j,k,l 

(33) 

^Kl(ti)nj^(tj) 

^I^kl^^i^^mn 

'V| 

= 0 (kl^^mn) 

(34) 

i.e.,  and  are  assumed  to  be  spatially  uncorrelated  as  well. 

Since  one  complete  scan  is  a measurement  array,  there  are  effectively 
64  measurements  available  at  each  sample  time. 
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III.  Extended  Kalman  Filter  Mcxlel 
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I 

i 

I 

1 
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Introduction 

I 

The  approach  used  in  the  design  of  the  Kalman  filter  was  to  keep 
the  models  for  target  motion  and  atmospheric  jitter  as  simple  as  pos- 
sible while  still  accounting  for  a time  correlation  in  the  dynamics  and 
a bandwidth  characteristic  of  the  jitter.  This  led  to  using  first  order 
Gauss-Markov  processes  along  each  independent  axis  direction,  to  describe 
both  the  target  dynamics  and  the  atmospheric  turbulence.  This  choice 
results  in  a four  state  filter  model  which,  with  the  proper  choice  of 
noise  strengths  and  correlation  times , will  be  descriptive  of  a wide 
rcuige  of  targets . 

State  Space  Model 

The  general  form  of  a stationary,  first  order,  Gauss-Markov  process 
can  be  expressed  as 


x(t)  = -1/t  x(t)  + w{t) 


(35) 


with 

E|w(t)j  = 0 

E|w{t)w(s)|  = ^ 6(t-s) 


(36) 

(37) 


where 

T = correlation  time 
w(t)  » white  noise  process 

a “ RMS  value  of  the  output  x 

Expressing  the  four  filter  states  (x^^,  x^^,  y^^,  y^)  in  this  manner  and 
putting  them  into  the  general  form  of 
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.aatBife. 


r 


^(t)  = Fp3^(t)  + Wj,(t) 


(38) 


where 


F 

i-F 


the  filter  plant  matrix 
3^(t)  = the  filter  state  vector 


w„(t)  = the  input  white  noise  vector 


gives  the  following  equation: 


.(t)  = 


■"/t. 


0 

0 

0 


-1/n 


0 

0 


0 

0 


0 

0 


0 

-V, 


3^{t)  + ^(t) 


(39) 


where 


E[^(t)]  = 0 


(40) 


and 


r 


i:|^^(t)^(s)J  = gp6(t-s) 


0 

0 

0- 


0 

^2 

0 

0 


0 

0 

0 


0 ^ 
0 
0 


^2  -» 


6(t-s) 


(41) 


State  Propagation 

The  filter  state  equations  shown  above  are  linear,  therefore,  stan- 
dard propagation  equations  can  be  used  to  propagate  the  filter  states 
between  measurement  times.  These  equations  are 


i^^i+1^ 


ip(ti^l,ti)i(tj) 


(42) 
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/^F(ti^l»T)Qp(T)$^(ti+i,T)  dX 


where  the  state  transition  matrix  ,t.)  satisfies  the  following 

~r  X 

differential  equation 


^(t  ,t^)  = Fp|p(t,ti)  (44) 

with  the  boundary  condition 

$„(t.,t.)  = I (the  identity  matrix)  (45) 

The  - eind  signs  refer  to  before  and  after  a measurement  is  tciken 
respectively . 

In  this  case,  the  plantmatrix  Fj,  is  a time-invariant  matrix, 
therefore,  the  ^(t,t^)  matrix  also  becomes  a constant  for  a given 
sampling  period. 


* , X F(t-ti) 

^(t,t^)  = e—  ^ 


-(t-tn) 
e 

Td 


(t-t.) 


- (t-t . ) 


-(t-t.) 
e — T*’ 


The  ^ matrix  is  also  a constant  matrix  for  all  time,  and  in  order  to 
get  the  proper  PMS  values  on  the  corresponding  processes  (Ref  4:193) , 


it  becomes 


0 


21 


0 


where 


variance  on  the  output  of  ttie  target 
dynamics  model 

variance  on  the  output  of  the  atmoshperic 
jitter  model 


The  integral  term  in  equation  (43)  now  becomes 


a^(l-e  ) 0 0 

D D 

0 0^  (l-e'^^VtA)  0 

A 

-2At/  , 

0 0 


2/-1  ”2At  / . 


(48) 


0 


0 


p(tj) 


r^-4t/TD 


-At/x 
e A 


e D 


l-e-2'^VT, 


a2(l-e-^“^/tA) 


(50) 


These  axe  the  equations  which  will  propagate  the  filter  states  between 
measurements , 

Measurements 

The  measurement  history  used  by  the  filter  to  help  predict  the 
pointing  directions,  will  be  the  outputs  from  the  forward  looking  infrared 
sensor  (FLIR) . These  outputs,  as  described  in  the  first  chapter,  are 
average  intensity  values  as  shown  in  Figure  1.  The  output  of  the  k-l-th 
pixel  can  be  described  mathematically  as 


^ (y-ypi^^)"]j<3^dy  + v,^3^(t.)  (51) 


where 


- the  actual  output  of  the  k-l-th  pixel  at  time  t^ 

= the  area  of  one  pixel 

= the  maximum  intensity  received  from  the  target 

= the  dispersion  of  the  Gaussian  intensity  function 
as  assumed  by  the  filter 

= coordinates  of  any  point  in  the  pixel  (variables 
of  integration) 


‘^oeak^V 


y e^(ti))  = the  center  of  the  intensity  pattern  with  respect 


to  the  center  of  the  FLIR  field  of  view  at  time  t^ 
additive  noise  term 


The  additive  noise  term  is  assumed  to  be  a spatially  and  temporarily 
white  Gaussian  noise  with  the  following  statistics : 


^O'kl^ti^]  = 0 

E[vj^l(ti)Vkl(tj)] 


(for  t^=tj) 
0 (for 


^tmn^^i^^kl^^j^]  = ° (for  all  ti,tj) 


(mn)j<(kl) 


It  is  used  to  account  for  euiy  uncertziinties  in  the  system,  such  as  back- 
ground noise , FLIR  noises , and  integration  approximations , 

In  general  state  space  form,  equation  (51)  is  of  the  form 


Zki(ti)  » h,^l(2L(ti).ti)  + v,,l(t^) 


1 


where  is  a nonlinear  function  of  the  states  at  time 

t^  . Because  of  the  8 by  8 tracking  array,  there  will  be  64  indepen- 
dent measurements  available  at  each  discrete  measurement  time.  These 
measurements  can  be  grouped  into  a 64-dimensional  vector  of  the  form 


where 


£(ti)  = h(x(tj^)  ,t^)  + v(tj^) 


E|v(ti)l  " - 

= Rpl 


(for  t^=hj) 
(for  tiJ^tj) 


Not  all  of  these  measurements  need  to  be  used  in  the  estimation  for  the 
peak  intensity.  In  fact,  a smaller  number  of  measurements  could  signifi- 
cantly reduce  computer  time  while  still  maintaining  good  tracking  ac- 
curacy . The  number  of  measurements  used  would  depend  on  two  factors . 
First,  it  would  depend  on  the  dispersion  of  the  Gaussian  intensity  func- 
tion O-  . For  a small  dispersion,  on  the  order  of  one  pixel,  most  of 
yp 

the  useful  information  is  contained  in  an  area  of  about  5 pixels  square . 
A smaller  dispersion  would  require  fewer  pixels  to  get  the  same  amount 
of  information.  The  second  factor  is  the  tracking  accuracy  needed. 

For  most  laser  applications,  high  accuracy  tracking  is  desired.  The 
best  accuracy  would  be  obtained  when  all  the  available  information  is 
used  in  the  prediction  process.  The  smaller  the  anount  of  information, 
the  Jess  accurate  the  prediction,  even  with  very  low  signal-to-noise 
ratios.  A third  factor  that  may  be  considered,  is  the  trade-off  be- 
tween estimation  performance  and  conputer  time  and  storage. 
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Measurement  Update  Equations 


In  order  to  handle  the  nonlinearities  in  the  measurement  equation, 
a special  form  of  the  Kalman  filter,  called  the  extended  Kalman  filter, 
was  used.  This  form  of  the  Kalman  filter  accounts  for  the  nonZineari- 
ties  by  linearizing  the  measurement  equation  about  the  most  recent  esti- 
mate of  the  states . This  linearization  process  is  reaccomplished  prior 
to  every  update,  thereby  maintaining  a current  set  of  equations  for 
updating  the  states . 

The  usual  form  of  the  measurement  update  equations  for  the  ex- 
tended Kalman  filter  (Ref  3:233)  are 


K(ti) 

(59) 

£(tt) 

= x(t^)  + K(tj^)  ,t^)  1 

(60) 

P(tt) 

= F(tT)  - K(t^)H(ti)P(tp 

(61) 

H(t^) 

_ 9h(x,t) 

^21 

Xf£(ti) 

(62) 

and  r(tj^)  is  the  64  dimensional  vector  of  actual  measurements. 

Notice,  however,  that  calculating  the  Kalman  filter  gains, 
k(tj^)  , requires  the  inversion  of  a 64  by  64  matrix  every  update  time. 
This  would  require  a large  amount  of  computer  time,  and  would  probably 
be  unfeasible  for  most  on-line  applications . Another  form  of  Ihe  up- 
date equations,  known  as  the  inverse  covarience  form  (Ref  3:257), 
eliminates  this  problem.  Using  this  form,  the  equations  become 

P“'(t^)  “ P~'(tJ)  + (63) 


\ 
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p(t+) 

- |£.-(4)|-‘ 

(64) 

K(ti) 

(65) 

x(tt) 

= x(t7)  + K(tj^)  {r(t^)-h(x(t^)  ,t~ ) } 

(66) 

where  the  matrix  and  the  r(tj^)  vector  are  defined  the  same  as 

above . Using  this  new  set  of  update  equations , only  two  4 by  4 inverses 
are  needed,  the  P^(tj^)  matrix,  and  the  P"Mtt)  matrix.  Normally, 
this  form  of  tlie  update  equation  would  also  involve  the  offline  inver- 
sion of  the  R(tj^)  matrix,  but  in  this  case,  the  pixel  noises  were 
assumed  independent  which  means  that 


R(t^)  = Rp£ 


(67) 


and  therefore. 


R^^(ti)  = (l/Rj.)I 


(68) 


Now  this  inversion  can  be  computed  once  offline  and  the  R~^(tj^) 
matrix  will  remain  constant  for  all  time. 

The  partial  derivative  matrix  shown  in  equation  (62)  can  be  easily 
computed  knowing  the  fact  that  the  integral  is  a linear  operator.  For 
this  reason , 


3h 


9 + I, 


9 


//^  "+^P-ypeak> 


20  2 

g 


3^x{tp 


Ap  92L 


exp{ 


-1 

20g2 


^P"ypeak) 


} 

dCdp 

3^£(ti) 


(69) 


(70) 
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where  the  partial  derivative  of  the  double  integral  term  is  equal  to 
the  double  integral  of  the  partial.  Expanding  this  into  matrix  form 


L 


and  evaluating  at  21  ~ results  in  a 64  by  4 matrix,  whose  ith 

row  is 

T 


j"  |(C“^eak^  *^f| 

//l(C-Xpeak^ ‘CpldCdp 
//l(p-?peak) ’CpldCdp 
//[(p-^peak^ -Cpld^dp J 


(71) 


for  x=(Xj^,  x^,  y^,  y^)  , , and  ?peak=yD+^A 


^ _ ‘^max 


Vg 


^^^2^  [^^“Vak)'  + (p-ypeak>']5  <^2) 
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Summary  of  Filter  Equations 

The  following  sets  of  equations  represent  the  extended  Kalman 
filter  algorithm  used  for  this  study . 

Propagation  Equations . 


where 


= 

^(ti^l)  x(t.) 

(73) 

= ; 

2D(t.) 

(74) 

e^^/'^D  0 

0 

0 

0 

0 

0 

•i+l'*^i^ 

0 0 

-At/T^ 
e D 

0 

0 0 

0 

(75) 

'mf 
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and 


Update  Equations . 


P“^(tt)  = PT'ttp  + H'^(tj|^)R"^(t^)H(t^)  (77) 

P{t^)  = [P“^(t]^)j“^  (78) 

K(ti)  = P(ti)H^{t^)R"Mti)  (79) 

x(tt)  = i (tT)  + K(t^) {r(t^)  - h(^(tT) ,t^) 1 (80) 


where  the  R ^ (tj^)  and  ^^_(tj^)  matrices  are  defined  as  in  equations 
(68)  cind  (69)  respectively. 

Now  that  the  truth  model  and  extended  Kalman  filter  model  have 
been  described,  the  next  chapter  will  outline  the  correlation  algorithm 
to  be  used  for  comparison  purposes , 
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IV . Correlation  Algorithm 


IntroductJ  on 

In  order  to  judge  the  effectiveness  of  the  extenled  Kalman  filter 
algorithm  in  tracking  a target,  its  performance  was  compared  to  that 
of  a correlation  algorithm  presently  being  tested.  Although  two  types 
of  correlation  algorithms  are  currently  being  examined  for  use  in 
operational  software  (Refs  4 and  5) , one  algorithm  was  readily  avail- 
able in  FORTRAN  source  code  (Ref  6) , and  has  been  tested  against  Hawk 
missile  test  flight  data.  This  algorithm,  therefore,  was  selected  as 
the  benchmark  for  coiaparison. 

This  chapter  focuses  on  the  method  used  by  the  selected  algorithm 
to  calculate  the  position  errors  needed  for  tracking.  A more  detailed 
analysis  can  be  found  in  Reference  4. 


Purpose 

The  purpose  of  any  tracking  algorithm  is  to  provide  estimates  cf 
the  position  of  the  target  with  respect  to  a reference  point,  usually, 
the  line  of  sight  of  the  sensor  tracking  the  target.  Knowing  this  in- 
fonration,  the  sensor  can  then  keep  the  target  centered  in  its  field 
of  view  by  appropriate  commands  to  a positioning  controller.  The  cor- 
relation tracker  estimates  these  relative  position  offsets  by  using 
correlation  coefficients  derived  from  the  cross-correlation  of  a target 
data  set  {t^},  and  a reference  data  set  {r^}.  This  correlation  co- 
efficient is  defined  as 


T TiRi 

i«»l 


Z Ti*  Z Ri* 

i=l  i=l 


(81) 
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i 


1 
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where  n is  the  number  of  sampled  data  points  in  each  of  the  data 
sets.  This  correlation  coefficient  is  a magnitude-independent  measure 
of  the  degree  of  resemblance  between  the  data  sets . This  means  that 
the  correlation  coefficient  is  a function  of  the  shape  of  the  two  data 
sets  and  is  independent  of  their  actual  magnitudes . 

The  correlation  technique  does  not  e^qsloit  any  knowledge  of  target 
dynamics,  nor  does  it  use  any  knowledge  of  disturbances  that  would 
cause  apparent  offsets  of  the  target,  so  that  these  effects  could  be 
filtered  out. 

One-Dimensional  Correlation 

The  following  discussion  will  be  a brief  overview  of  correlation 
in  one  dimension  as  presented  in  an  Aeronutronic  Ford  Corporation  re- 
port (Ref  4) . Let  {r^}  be  a set  of  uniformly  spaced  samples  of  a con- 
tinuous time  radiation  intensity  function,  received  from  a one-dimen- 
sional scan  of  a FLIR  sensor  across  a point  target.  Let  {t^}  be  a 
new  sample  set  taken  a short  time  later,  where  the  intensity  data  is 
assumed  to  have  the  same  shape,  but  may  have  been  translated  relative 
to  the  first  data  set  because  of  target  motion  in  the  dimension  of 
interest.  The  {r^}  and  {Tj^}  data  sets  can  be  expressed  as 

{r^}  = sampled  [s(x)]  (82) 

{Tj^}  = sampled  [s(x  + Ax)]  (83) 

where  S(x)  is  the  continuous  time  intensity  function  and  Ax  is  the 
intensity  offset  on  the  FLIR  image  plane.  E:q>anding  S(x  + Ax)  in 
a taylor  series  gives 

S(x  + Ax)  = S(x)  + 8s(x)  A X + H.O.T.  (84) 

3x 
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Now  by  dropping  the  higher  order  terms,  {t^}  can  be  approximated  as 


This  equation  represents  one  equation  in  three  unknowns  (P,  Ax,  sign 
of  Ax) . Two  more  equations  can  be  established  by  shifting  one  of  the 
sampled  data  sets  left  and  right  by  one  sample.  This  results  in  the 
following  three  equations : 


S' 

= 1 

- P^Ax^ 

(92) 

s' 

= 1 

- P^(l  - Ax)^ 

(93) 

s' 

= 1 

- P^(l  + Ax)^ 

(94) 

where  Cp  represents  the  correlation  coefficient  produced  when  the 
two  data  sets  are  aligned  to  the  closest  sample . and  C3  repre- 

sent the  correlation  coefficients  obtained  when  the  data  set  is  shifted 
left  and  right  by  one  sample  respectively.  Solving  these  three  equations 
simultaneously  gives 


Ax 


(95) 


This  equation  represents  the  offset,  to  a fraction  of  a sample,  between 
the  data  sets,  and  is  a direct  measure  of  the  positional  change  of  the 
target  over  the  sampling  time. 

The  actual  equation  for  Ax  used  for  this  study  was 

Ax  = l/Cp^  - 1/03^ 

2 < ^ + 1 - 2 \ 

cp  cp"  c^2  (96) 


This  equation  was  developed  at  the  Air  Force  Weapons  Laboratory  after 
the  equations  leading  up  to  equation  (95)  were  rederived,  eliminating 
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the  simplifying  assumption 


1 1 + P^Ax^  - P^Ax^  ~ 1 _ tj2Av.2 

1 + P^Ax2  ■ 1 + p2Ax2  (97) 

The  derivation  of  equation  (96)  can  be  found  in  a memo  from  the  Air 
Force  Weapons  Laboratory  to  the  Ford  Corporation  (Ref  7) . 

In  order  to  present  a clearer  representation  of  what  is  taking 
place , let  Figure  6 represent  a one-dimensional  Gaussian  intensity 
profile  of  a point  source  target  moving  horizontally  to  the  right. 

At  time  tj  , the  target  is  centered  in  the  field  of  view  of  eight 
horizontal  pixels,  and  at  time  t2  , the  target  has  moved  some  dis- 
tance r . Assuming  ideal  conditions  (no  change  of  shape  and  perfect 
measurement  data) , the  intensity  profile  will  also  move  to  the  right 
some  distance  d . Figures  (7a)  and  (7b)  represent  the  discretized 
value  of  the  intensity  profile  at  time  t2^  and  t2  • Each  discretized 
value  represents  the  averaged  sampled  output  of  the  detectors  as  they 
cire  scanned  across  the  target.  These  discrete  values  become  the  ref- 
erence and  target  data  sets  |Ri)  and  {Tj^l  respectively.  The  center  of 
the  intensity  function,  to  the  nearest  pixel  in  each  data  set,  is  esti- 
mated by  the  simple  mean  approximation. 


J(t) 


8 

E Ii(t)*i 
i=l 


INT 


8 


E Ii(t) 
i=l 


(98) 


where 

Jit)  = ein  integer  from  one  to  eight  representing  the  pixel , 
from  left  to  right,  closest  to  the  center  of  the 
intensity  function 


2 3 4 5 6 7 8 
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Pixels  Pixels 

Discrete  Representation  at  tj^  b.  Discrete  Representation  at  t. 

Figure  7.  Discrete  Representations  of  Intensity  Profile 


and 


= the  intensity  value  in  the  ith  pixel  at  time  t ; 


INT(*)  = a function  which  rounds  the  argviment  to  the  nearest 
integer. 

The  offset  of  the  intensity  profile  to  the  nearest  pixel  can  now  be  esti- 
mated as  the  quantity  J(t2)  - JCt^^)  . The  reference  array  {r^}  , which 
is  the  array  of  values  at  time  tj^  , is  now  shifted  by  J(t2  - J{t2) 
pixels.  Those  values  tliat  are  shifted  out  of  the  array  are  dropped,  and 
the  values  at  the  opposite  end  of  the  array  are  replaced  by  zeros . An 
example  of  this  shifting  process  is  shown  in  Figure  8.  The  three  corre- 
lation coefficients  > ^2  • ^3  • ^ computed  for  the 

shifted  reference  array  {Ri}s  target  array  {t^}.  Equation  (96) 

is  then  used  to  estimate  the  remaining  offset  to  a fraction  of  one  pixel . 
The  total  offset,  d , is  then  estimated  by  the  equation 

d = J(t2)  - J(ti)  + Ax  (99) 


Two-Dimensional  Correlation 

For  a two-dimensional  intensity  function,  the  above  process  is 
performed  for  every  row  and  then  every  colvunn  in  the  8x8  arrays. 

This  will  result  in  eight,  generally  different,  horizontal  offsets,  and 
eight  vertical  offsets.  Because  of  this  structure,  the  image  position 
offsets  in  the  horizontal  and  vertical  directions  are  arrived  at  inde- 
pendently. The  row  correlation  determines  the  horizontal  offset,  while 
the  column  correlation  determines  the  vertical  offset.  If  the  maximum 
correlation  coefficient  C2  for  each  row  or  column  is  less  than  some 
preset  value,  that  row  or  column  is  considered  as  containing  poor 
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{Rj^}  = 12464211 


{r.  }„  = 0 0 1 2 4 6 4 2 

i s 


a.  Reference  Array  Before  Shift 


b.  Reference  Array  After  Shifting 
Two  Pixels  Right 


Figure  8.  Example  of  an  Array  Shift 


information,  and  the  computed  offset  is  discarded.  The  remaining  row 
and  column  offset  values  are  then  averaged  separately.  The  average 
of  the  offset  values  for  the  rows  is  the  overall  estimate  of  target 
position  offset  in  the  horizontal  direction,  and  the  average  of  the 
offset  values  for  the  columns  is  the  estimated  target  position  offset 
in  the  vertical  direction. 

This  simple  averaging  technique  is  a simplification  of  the  method 
used  by  the  Rod  Corporation  in  their  report  (Ref  8) . The  following 
discussion  explains  the  reason  for  this  simplification.  The  Aero- 
nutronic  Ford  Corporation  algorithm  used  a least  squares  line  fit  of 
the  values  in  each  direction.  One  line  was  fit  to  the  individual  row 
offsets,  and  one  to  the  column  offsets.  The  total  image  offset  in 
each  direction  was  then  computed  in  a way  unspecified  in  the  available 
sources.  One  method  would  be  to  determine  the  intersection  of  each 
line  with  the  appropriate  x or  y axis . Another  might  be  to  deter- 
mine the  intersection  of  the  two  regression  lines . The  slope  of  these 
lines  is  then  a measure  of  the  target  rotation  over  that  sample  period. 

In  this  study,  however,  the  GaussicUi  intensity  function  is  not 
directly  affected  by  target  rotation.  For  this  reason,  the  slope 
will  contain  no  useful  information  and  would,  in  fact,  be  either  zero 
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or  ninety  degrees  under  ideal  conditions , where  the  offsets  in  each 
direction  will  all  be  equal.  This  entire  method  was,  therefore, 
simplified  to  simply  finding  the  center  of  mass  of  the  offset  points 
in  each  direction,  defined  as  tlie  average  point  since  all  points  are 
equally  weighted. 

Summary 

The  correlation  algorithm  described  above,  is  only  one  of  several 
algorithms  that  could  be  used  in  a correlation  tracker.  This  algorithm 
was  readily  availabJe,  and  was  suggested  for  use  by  the  Air  Fbrce  Weapons 
Laboratory.  No  attempt  was  made  to  determine  if  it  was  the  most  effi- 
cient or  optimal  in  any  sense.  Care  was  taken  only  to  make  certain 
that  both  this  algorithm,  and  the  extended  Kalman  filter  algorithm  re- 
ceived the  same  inf onoation  from  the  PLIR  sensor  model , In  this  way , 
the  filter  algorithm  could  be  compared  for  accuracy,  in  a statistical 
sense , against  an  algorithm  already  in  use , to  examine  the  utility  of 
an  extended  Kalman  filter  in  the  tracking  loop. 
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V.  Computer  Simulation 

, > 

‘ \ 

Introduction 

The  computer  used  to  simulate  the  models  of  the  previous  three 
chapters  was  the  CDC-6600.  As  with  any  computer  simulation,  the  time 
and  memory  requirements  vary  directly  with  the  complexity  of  the 
problem.  The  key,  of  course,  is  to  minimize  these  requirements  and 
still  achieve  an  algorithm  that  will  perform  well  in  a real  world 
environment.  This  chapter  will  discuss  the  methods  used  to  arrive 
at  this  balance . 

Program  Modularity  and  Efficiency 

An  attempt  was  made  to  keep  the  program  as  modular  as  possible 
by  using  subroutines  for  major  calculations.  The  purpose  of  the  study 
was  not,  however,  to  create  a generalized  program  to  Work  under  all 
possible  variations  of  filter  design.  For  this  reason,  changing  any 
of  the  models  discussed  previously  would  probably  necessitate  at  least 
some  modification  to  the  program  modules.  The  FORTRAN  source  code  of 
the  program  used,  is  liste^d  in  Appendix  D. 

The  large  word  size  eind  relative  speed  of  the  CDC-6600,  compared 
to  a machine  used  for  on-line  implementation,  also  compensates  for 
many  inefficiencies  in  programming  techniques.  Although  the  algorithms 
developed  earlier  are  meant  to  be  implemented  on-line  eventually,  the 
program  used  would  have  to  be  restructured  to  make  it  as  efficient  and 
numerically  precise  as  possible . 

Approximations 

One  of  the  most  important  aspects  of  the  computer  simulation  was 
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determining  a method  for  evaluating  the  outputs  of  the  FLIR  sensor. 
Mathematically,  this  would  involve  evaluating  the  expression 

T"  ^max  \ - Vak^^  + (y  - Ypeak)  ^dxdy  (100) 

*p  i-j-th  2ag^  1 ^ x-  j 

pixel 

as  described  in  equation  (30)  for  the  truth  model.  Similar  expressions 
would  have  to  be  evaluated  for  the  filter  model  to  arrive  at  the  fil- 
ter's estimate  of  the  measurements,  equation  (51),  and  to  evaluate  the 
partial  matrix  found  in  equation  (71) . 

E^guation  (100)  represents  the  average  value  of  the  two-dimensional 
Gaussian  intensity  function 

W 1 - "^eak)  ^ + (y  - ypeak>  ^ | ^ 

within  the  pixel  of  interest.  A reasonable  approximation,  for  on-line 
usage,  to  equation  (100)  would  be  siinply  to  evaluate  the  intensity 
function  at  the  midpoint  of  the  pixel  of  interest.  A better  approxi- 
mation would  be  to  subdivide  the  pixel  into  smaller  areas , perform  the 
s^lIne  midpoint  evaluations  for  each  of  these  smaller  areas , and  to 
average  the  result,  obviously,  the  more  subdividions  that  are  made, 
the  better  the  approximation. 

Because  the  truth  model  is  designed  to  be  a good  representation 
of  the  real  world,  an  accurate  evaluation  of  the  integral  is  needed. 

To  meet  this  requirement,  each  pixel  was  split  into  sixteen  equal  sub- 
divisions. Tlie  filter,  on  the  other  hand,  used  the  intensity  level 
at  a single  midpoint  to  approximate  the  integral. 

White  Noise  Terms 

Throughout  the  simulation,  there  was  a need  to  produce  discrete 
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values  of  a Gaussian  distributed,  white  sequence  with  a given  variance. 


In  order  to  create  these  values , the  RAND  function  of  the  FORTRAN 
library  was  used.  Tliis  function  is  a random  number  generator  with  a 
uniform  probability  density  function  between  0 and  1.  Twelve  inde- 
pendent random  samples  from  this  function  were  processed  with  tlie 
algorithm 

12 

v{k)  = a„{  I ull2(k  - 1)  + i]  - 6}  (102) 

i=l 

for  k = 1 , 2 , 3 , ... 
where 

v(k)  = a discrete  realization  of  the  white,  Gaussian 
sequence ; 

= the  standard  deviation  of  the  desired  white  noise; 

u(k)  = random  nvimber  from  the  RAllD  function; 
to  generate  a single  sample  of  the  desired  white,  Gaussian  sequence. 
The  above  algorithm  produces  a sequence  of  numbers  that  is  very  nearly 
Gaussian  distributed  with  zero  mean  and  the  desired  variance.  This  is 
triie  because  of  the  concepts  set  forth  in  the  Central  Limit  Theorem 
(Ref  3:103)  . 

Monte  Carlo  Analysis 

Determining  how  any  computer  algorithm  will  perform  in  a real 
world  environment  is  a difficult  process.  Unlike  the  carefully  con- 
trolled computer  simulation,  the  actual  environment  is  constantly 
changing.  No  two  experiments  are  ever  the  same,  cind  to  test  all  pos- 
sible cases  would  be  an  impossible  task.  For  this  reason,  a Mante 
Carlo  analysis  is  used  to  accomplish  a performance  evaluation.  This 
method  actually  constructs  a seunple-by-scunple  simulation  using  random 
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number  generators  to  produce  the  outputs  of  random  noise  sources  in 
the  simulation.  Tlie  noisy  system  measurements  are  then  processed 
by  the  extended  Kalman  filter,  or  the  correlator  algorithm,  to  gen- 
erate state  estimates.  Many  of  these  simulations  are  performed  and 
the  Sdimple  statistics  of  the  error  between  the  state  estimates  and 
the  states  of  the  truth  model  are  then  computed.  The  trend  of  these 
error  statistics  shows  how  the  estimation  algorithm  performs  in  a 
general  sense.  Ideally,  a larger  number  of  simulations  should  be  per- 
formed to  make  sure  that  the  error  sample  statistics  are  indeed  accu- 
rate representations  of  the  true  statistical  characteristics . Real- 
istically, however,  computer  time  requirements  dictate  a smaller 
number.  One  .tiecins  to  estimate  the  number  of  runs  needed  is  to  plot 
the  standard  deviation  of  the  error  at  a number  of  chosen  times  versus 
the  lumber  of  simulations . As  the  number  of  simulations  increase , the 
stcindard  deviation  should  converge  to  a steady  value.  This  convergence 
shows  that  the  error  statistics  are  being  adequately  described.  An 
example  of  this  convergence  is  shown  in  Figure  9 . It  shows  the  vari- 
ance characteristics  of  a representative  error  at  four  different  times. 
It  also  shews  that  20  simulations  were  adequate  for  this  case.  These 
types  of  plots  were  produced  for  every  case  studied  to  make  sure  the 
statistics  converged.  Twenty  simulations  were  found  to  be  adequate 
in  all  cases . 

Tuning  the  Extended  Kalman  Filter 

Tuning  the  filter  is  an  off-line  process  to  determine  the  proper 
noise  strengths  and  correlation  times  so  the  filter  will  perform  well 
in  any  situation.  This  is  necessary  because  the  filter  is  purposely 
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designed  to  be  a reduced  order,  simplified  model  of  the  real  world, 
and  must  overcome  the  effects  of  non-whiteness,  non-Gaussian  noises, 
nonlinearities,  and  approximations  to  a true  environment. 

The  tuning  process  is  done  by  varying  the  values  of  the  strengths 
of  the  white  noise  terras  in  both  the  differential  state  equation  and 
the  filter  measurement  equation.  In  this  study,  this  corresponds  to 
altering  the  values  of  the  Q^,  matrix  given  in  equation  (47)  and  the 
value  of  Rp  in  equation  (50).  Cl'ianging  these  values,  along  with  the 
initial  conditions,  will  cause  the  statistics  of  the  error  committed 
by  the  filter  in  estimating  the  states  to  vary.  The  objective  of 
tuning  is  to  vary  these  strengths  in  order  to  minimize  the  estimation 
error  in  some  appropriate  manner.  Here,  the  object  is  to  minimize 
the  error  on  all  filter  states  simultaneously.  This  minimization  is 
totally  problem  dependent,  however,  and  other  problems  may  call  for 
minimizing  only  a given  number  of  states,  using  a weighted  average, 
or  using  some  type  of  cost  function  to  measure  the  minimum. 

One  method  used , in  many  cases , to  aid  the  tuning  process , is 
to  compare  the  filter's  estimate  of  the  variance  of  the  error  with 
the  actual  error  variance  for  each  state.  When  these  values  match 
well  over  the  time  interval  examined,  then  the  actual  error  committed 
by  the  filter  should  be  minimized.  If  the  filter  is  underestimating 
or  overestimating  the  true  variance,  the  filter  geiins  will  not  be 
optimal,  and  will  result  in  a larger  true  error.  It  is  important, 
therefore,  that  the  truth  model  be  as  accurate  as  possible.  The  fil- 
ter will  work  well  on-line  only  if  the  truth  model  is  an  accurate 
representation  of  the  real  world. 
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VI . Results 

Introduction 

This  chapter  presents  the  results  of  several  computer  simulations 
designed  to  evaluate  the  performance  of  the  extended  Kalman  filter 
algorithm  in  the  open-loop  tracking  problem.  As  requested  by  the  Air 
Force  Weapons  Laboratory,  the  performance  of  the  extended  Kalman  filter 
was  compared  to  the  performance  of  the  correlation  tracker  in  five 
specific  areas : 

- performance  as  a function  of  signal-to-noise  ratio  where 

S ^max  ~ nisan  background  noise 
N RMS  value  of  background  noise 

- performance  as  a function  of  intensity  pattern  size 
(beam  width) ; 

- performance  as  a function  of  the  ratio  of  RMS  target 
motion  to  RMS  atmospheric  jitter; 

- performance  as  a function  of  target  correlation  time; 

- performance  as  a function  of  mismatches  between  the  true 
intensity  pattern  and  that  assumed  by  the  filter. 

Due  to  the  limitations , all  of  the  above  areas  were  covered  with  the 

first  three  areas  receiving  the  most  detailed  analysis . 

Correlation  Algorithm 

The  correlation  tracking  algorithm  (Ref  6)  was  modified  slightly 
after  several  problems  were  discovered  during  the  study . The  first 
problem  occurred  in  the  LOCATE  subroutine  (Appendix  D)  when  a line  of 
reference  data  was  shifted  by  more  than  8 pixels , resulting  in  an  entire 
line  of  zeros  and  an  indefinite  correlation  coefficient.  To  solve  this 
problem,  a check  was  put  into  the  subroutine  to  assure  that  each  row 
or  column  was  limited  to  a shift  of  less  than  8 pixels  when  trying  to 
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align  the  target  and  reference  data  sets.  If  the  target  and  reference 
rows  or  columns  were  not  aligned  after  7 pixel  shifts  in  the  same 
direction,  then  the  data  from  that  row  or  column  was  considered  "poor" 
and  was  discarded. 

The  second  modification  involved  a control  in  the  DRIVER  program 
designed  to  reject  low  quality  offset  data  by  checking  the  maximum 
correlation  coefficient  for  each  row  and  column  against  a preset  cut- 
off value.  In  the  routine  obtained  from  the  Weapons  Laboratory,  the 
row  correlation  loop  eliminated  offset  data  with  a maximum  correlation 
coefficient  less  than  0.95,  while  the  column  loop  used  0.90  as  the 
cutoff  value.  There  was  no  documented  reason  for  this  difference,  nor 
was  there  any  justification  for  these  seemingly  arbitrary  cutoff  values. 
It  was  discovered,  however,  that  changing  this  cutoff  value  led  to 
drastic  differences  in  the  performance  of  the  correlation  algorithm. 
These  differences  are  documented  and  discussed  later  in  this  chapter. 

For  these  reasons,  this  study  used  0.95  as  the  cutoff  for  both  row  and 
column  correlation.  Similar  controls  in  the  LOCATE  subroutine  were 
eliminated  since  they  were  redundant. 

Cases  Studied 

In  order  to  investigate  the  five  specified  areas  in  the  time  avail- 
able, the  sixteen  scenarios  listed  in  Table  1 were  simulated.  As  in- 
structed by  the  Air  Force  Weapons  Laboratory,  the  specific  values  used 
in  each  simulation  were  combinations  of  the  following; 

- signal-to-noise  (S/N)  ratios  of  20,  10,  1 

- RiMS  target/atmospheric  motion  ratios  (o^/a^)  of  5,  1,  .2 

- Gaussicui  bean  dispersions  (a_  ) of  3,  1 pixels 

gip 
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Table  I 

List  of  Cases  studied 


Case 

^T 

S/N 

°gF 

1 

1. 

20. 

1.0 

0.2 

5.0 

3. 

3. 

2 

1, 

20. 

1.0 

1.0 

1.0 

3. 

3. 

3 

1. 

20. 

0.2 

1.0 

0.2 

3. 

3. 

4 

1. 

20. 

1.0 

0.2 

5.0 

1. 

1. 

5 

1. 

10. 

1.0 

0.2 

5.0 

3. 

3. 

6 

1. 

10. 

1.0 

1.0 

1.0 

3. 

3. 

7 

1. 

10. 

0.2 

1.0 

O..** 

3. 

3. 

8 

1. 

10. 

1.0 

0.2 

5.0 

1. 

1. 

9 

1. 

1. 

1.0 

0.2 

5.0 

3. 

3. 

10 

1. 

1. 

1.0 

1.0 

1.0 

3. 

3. 

11 

1. 

1. 

0.2 

1.0 

0.2 

3. 

3. 

12 

1. 

1. 

1.0 

0.2 

5.0 

1. 

1. 

13 

5. 

20. 

2.0 

0.4 

5.0 

3. 

3. 

14 

1. 

20. 

2.0 

0.4 

5.0 

3. 

3. 

15 

1. 

20. 

1.0 

0.2 

5.0 

3. 

1. 

16 

1. 

20. 

3.0 

0.6 

5.0 

1. 

1. 

I 
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- Target  correlation  times  of  1,  5 seconds 

The  value  specified  by  a_  in  the  table  is  the  value  for  the  Gaussian 
dispersion  assumed  by  the  filter. 

The  sixteen  cases  shown  in  Table  I were  chosen  to  allow  a general 
framework  in  which  direct  comparisons  could  be  made.  Cases  1 through 
12  compare  four  basic  scenarios  using  the  three  signal-to-noise  ratios 
and  the  two  values  of  beam  dispersion  shown  above.  These  four  basic 
scenarios  allow  direct  comparison  of  the  performance  of  the  Kalman 
filter  under  various  atmospheric  conditions . The  remaining  four  cases 
were  used  to  show  the  effects  of  using  a large  value  for  target  correla- 
tion time,  determining  the  effects  of  filter  mismatch,  and  determining 
the  effects  of  using  a large  RMS  target  motion  combined  with  a small 
target  pattern  size.  The  results  of  the  simulations  performed  are  pre- 
sented in  Appendix  E as  mean  ±lo  error  plots  for  the  correlator,  the 
filter  estimate  of  target  dynamics,  and  the  filter  estimate  of  atmo- 
spheric jitter.  The  mean  and  sigma  values  are  sample  statistics  gen- 
erated from  a Monte  Carlo  simulation  using  20  runs . Because  of  the 
similarity  between  the  horizontal  (x)  and  vertical  (y)  channels,  only 
the  mean  ±la  error  plots  for  the  horizontal  channel  are  presented. 

Tuning 

Because  of  the  similarity  of  the  truth  and  filter  models  used  in 
this  study,  timing  of  the  extended  Kalman  filter  was  relatively  easy. 

The  only  difference  between  the  models  was  the  set  of  equations  used 
to  simulate  atmospheric  disturbance . whereas  the  truth  model  used  a 
third  order  Markov  process  to  simulate  the  atmospheric  jitter,  tlie 
extended  Kalman  filter  used  a first  order  process  with  a break  frequency 


48 


) 

t 

I 

» 

set  equal  to  the  lower  breeik  frequency  of  the  higher  order  model. 

Similarly,  the  varieince  for  the  white  noise  process  used  in  the  filter 
measurement  equation  (53)  was  set  equal  to  the  sum  of  the  variances 
of  the  processes  simulating  the  PLIR  and  background  noises  in  equa- 
tions (32)  and  (33) . This  is  appropriate  since  all  of  the  above  pro- 
cesses are  simulated  as  independent,  zero  mean,  and  Gaussian.  Since 
the  value  of  Rp  is  now  fixed,  adequate  tuning  results  when  the  strengths 
of  the  white  noise  terms  and  the  correlation  times  for  both  the  truth 
and  filter  models  for  dynamics  are  set  equal , and  when  the  RMS  values 
for  atmospheric  jitter  for  both  models  are  set  equal.  For  the  third 
order  atmospheric  model,  this  mecins  that  the  correct  value  of  the  gain, 

K,  must  be  calculated  to  achieve  the  desired  RMS  value  on  the  output 
process . 

Figures  10,  11,  and  12  show  plots  of  the  actual  versus  the  filter- 
computed  value  of  the  RMS  error  committed  in  estimating  target  motion 
in  the  horizontal  direction  for  cases  2,  6,  and  10  after  the  above 
tuning  was  performed.  Figures  13,  14,  and  15  show  the  Scime  type  of 
plots  for  the  atmospheric  jitter.  The  smoother  curve  represents  the 
filter  estimate  of  the  RMS  error.  The  plots  show  typical  results  for 
signal-to-noise  ratios  of  20,  10,  and  1 respectively.  As  shewn  in 
Figures  12  and  15 , the  filter  begins  to  underestimate  the  RMS  error 
as  the  signal-to-noise  ratio  gets  low.  This  is  due  to  the  effects  of 
the  small  mismatch  between  the  truth  model  and  the  filter  finally 
becoming  apparent.  Because  of  the  way  the  gp  matrix  and  the  value 
for  Rp  are  chosen  in  the  above  tuning  process , the  tuning  was  independent 
of  the  signal-to-noise  ratio  and  the  Gaussian  beam  dispersion.  This 
allows  comparison  of  the  outputs  of  identically  tuned  filters  for 
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Figure  10.  FILTER  VS.  RCTUHL  SIGMR  PLOT  ( S/N  = 20) 
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different  signal-to-noise  ratios  and  different  beam  dispersions . 

These  comparisons  will  show  how  the  performance  of  the  filter  will 
change  as  the  tracking  environment  changes. 

Discussion 

In  comparing  the  performance  of  the  extended  Kalmein  filter  to 
the  correlation  algoritlim,  the  five  specific  areas  mentioned  earlier 
were  excunined.  The  extended  Kalman  filter  algorithm  performed  ex- 
tremely well  as  the  signal-to-noise  ratio  was  lowered,  with  no  notice- 
able change  between  the  ratios  of  20  and  10,  and  only  a slight  degra- 
dation in  performarice  shown  at  a signal-to-noise  ratio  of  1.  Tables 
II  and  III  present  a comparisc  n of  the  performance  of  the  two  algorithms 
in  mean  and  1 sigma  tracking  error  for  the  three  signal-to-noise  ratios 
examined.  The  numbers  in  the  tables  are  values  that  were  estimated 
from  the  applicable  plots  and  are  only  good  to  1 or  2 significant 
figures.  Table  II  represents  the  average  of  the  estimated  values  for 
cases  1 through  12  where  the  beam  dispersion  was  set  to  3 pixels, 
while  Table  III  represents  the  estimated  values  when  the  beeim  disper- 
sion was  set  to  1 pixel.  The  tebles  also  represent  values  estimated 
at  the  end  of  tlic  5 second  run  time  since  this  is  where  the  maximum 
errors  occurred  in  most  cases . As  shown  by  Uie  plots  for  cases  1 
through  12,  both  algorithms  started  out  with  low  errors  due  to  the 
chosen  initial  conditions  of  zero  position  error.  The  correlation 
algorithm  quickly  diverged  in  almost  all  cases,  resulting  in  signifi- 
cant errors  after  only  a few  sec'uds.  The  extended  Kalman  filter,  on 
the  other  hand,  had  mean  errors  of  close  to  zero  for  all  signal-to- 
noise  ratios,  and  although  the  1 sigma  error  value  rose  from  .2  to  .8 
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Table  II 

Mean  Error  and  la  Error  Comparisons  for  the  Correlation 


Tracker  and  the 

Extended  Kalman 

Filter  with  Og=3 

Pixels 

s 

Correlation 

Tracker 

Extended  Kalman  Filter 

N 

mean  error 
(pixels) 

10  error 
(pixels) 

mean  error 
(pixels) 

10  error 
(pixels) 

20 

0.5 

1.5 

0.0 

0.2 

10 

3.0 

3.0 

0.0 

0.2 

1 

15.0 

30.0 

0.0 

0.8 

Table  III 

Mean  Error  and  lo  Error  Comparisons  for  the  Correlation 
Tracker  and  the  Extended  Kalmari  Filter  with  dg=l  Pixel 


S 

N 

Correlation  Tracker 

Extended  Kalman  Filter 

mean  error 
(pixels) 

lo  error 
(pixels) 

mean  error 
(pixels) 

lo  error 
(pixels) 

20 

7.0 

8.0 

0.0 

0.2 

10 

8.0 

10.0 

0.0 

0.2 

1 

15.0 

30.0 

0.0 

0.8 

f 

t 
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pixels  as  the  signal-to-iioise  ratio  decreased,  this  difference  is  not 
significant,  when  the  dispersion  of  the  Gaussian  intensity  function 
was  decreased  to  1 pixel , the  performance  of  the  correlation  tracker 
was  significantly  degraded  with  errors  ranging  from  7 to  10  pixels. 

Tliis  same  change  had  essentially  no  effect  on  the  performance  of  the 
extended  Kalman  filter  algorithm.  This  difference  in  performance  was 
probably  due  to  the  a priori  knowledge  of  beam  width  and  target  motion 
assumed  by  tlie  Kalman  filter.  As  the  size  of  the  beam  width  is  de- 
creased, the  useful  tracking  information  is  spread  over  a smaller 
number  of  pixels . This  knowledge  can  be  exploited  by  the  filter  to 
give  better  performance,  especially  at  the  lov/er  signal-to-noise 
ratios.  This  suggests  that  using  data  from  a smaller  region  about 
the  projected  target  location  may  be  more  efficient  than  using  the 
entire  8 by  8 array.  Even  the  filter  performance,  however,  is  slightly 
degraded  at  a low  signal-to-noise  ratio  since  t}ie  approximations  used 
in  calculating  the  H matrix  and  in  performing  the  integration  are  not 
as  accurate . 

When  the  ratio  of  BMS  target  motion  to  RMS  atmospheric  jitter  was 
decreased  from  5 to  1 , the  filter  showed  a corresponding  increase  in 
actual  error  with  a mean  still  close  to  zero.  Comparing  cases  1 and 
2 shows  that  the  1 sigma  value  on  the  error  jumps  from  about  .2  pixels 
to  about  .5  pixels.  This  trend,  however,  is  reversed  as  the  ratio  is 
decreased  further  as  shown  by  the  output  of  case  4 , where  the  ratio 
was  equal  to  ,2.  Here,  tJie  1 sigma  value  for  the  error  goes  back  down 
to  about  .2  pixels  again.  This  seems  to  imply  that  the  filter  can 
distinguish  between  the  true  and  apparent  motion  easier  when  there  is 
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a large  difference  between  the  two.  This  type  of  trend  was  harder  to 
distinguish  in  the  outputs  of  the  correlation  algorithm  due  to  the 
divergent  characteristics  of  the  statistics . A comparison  of  the  same 
cases,  however,  seem  to  present  no  apparent  difference  in  accuracy  due 
to  a chcmge  in  ratio.  In  all  three  cases,  the  correlation  algorithm 
had  a 1 sigma  error  value  of  around  1.  pixel.  Even  in  the  worst  case, 
therefore,  the  extended  Kalmain  filter  algorithm  performed  better. 

Increasing  the  correlation  time  for  the  target  in  both  the  truth 
model  and  the  filter  from  1 to  5 seconds  had  no  discernable  effect  on 
the  performance  of  either  the  extended  Kalman  filter  algorithm  or  the 
correlation  algorithm.  A comparison  of  the  output  for  cases  13  and  14 
show  almost  no  chcuige  in  tlie  statistics.  Tliis  is  probably  due  to  the 
small  difference  in  correlation  tines  used.  If  the  correlation  time 

were  changed  to  a large  value,  on  th«-  ■ r of  1 minute  or  more,  the 

•*  • 

tracking  accuracy  should  improve  since  high  frequency  jitter  term 

could  be  filtered  more  easily  from  the  overall  motion  apparent  in  the 
target  plane.  A much  shorter  correlation  time,  on  the  other  haind, 
would  make  the  actual  target  motion  appear  more  like  a white  noise 
texrm  and  could  introduce  a larger  error. 

Case  15  examined  the  effects  of  mismatching  the  shape  of  the  in- 
tensity function  between  the  truth  model  and  the  filter.  The  filter 
assumed  a dispersion  on  the  Gaussian  intensity  function  of  1 pixel 
and  the  truth  model  value  was  set  at  3 pixels . This  scenario  was 
designed  to  test  the  sensitivity  of  the  filter  to  poor  assumptions  in 
target  shape.  In  this  case  the  filter  appeared  to  be  rather  robust, 
resulting  in  a mean  error  of  cibout  zero  with  a 1 sigma  value  cf  only 
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.3  pixels.  I’lhen  compared  to  case  1,  where  the  filter  assumed  the 
correct  value  of  3 pixels  for  the  dispersion,  there  was  only  a slight 
degradation  in  performance . Case  1 resulted  in  a mean  error  of  about 
0.  with  a 1 sigma  value  of  .2  pixels.  Although  results  based  on  one 
case  are  inconclusive , this  shows  that  the  filter  could  be  relatively 
insensitive  to  poor  assumptions  and  further  testing,  such  as  mis- 
matches in  target  shape,  target  motion,  or  signal-to-noise  ratios, 
would  be  warranted. 

The  worst  perforraance  of  the  extended  Kalman  filter  resulted 
when  the  value  of  RMS  target  motion  was  very  high , as  in  case  16  where 
it  was  set  to  3 pixels . This  resulted  in  RMS  tracking  errors  of  2 
to  3 pixels  after  5 seconds  of  run  time.  This  is  not  surprising,  how- 
ever, since  this  magnitude  of  target  motion  forced  the  target  entirely 
out  of  tlie  field  of  view  of  the  fLIR  on  many  passes.  This,  of  course, 
resulted  in  extremely  poor  measurement  data  which  would  tend  to  drive 
the  tracking  errors  up,  and  suggests  a larger  field  of  view  for  update. 

One  interesting  result  is  the  divergence  characteristic  of  the 
correlation  tracker.  This  divergence  was  noted  in  almost  all  cases, 
even  with  high  signal-to-noise  ratios.  It  was  discovered,  however, 
that  in  many  cases  this  divergence  could  be  reduced  or  even  eliminated 
by  increasing  the  cutoff  control  values  which  reject  the  offset  data 
from  poorly  correlated  rows  and  columns.  This  reduction  is  distinctly 
shown  in  Figures  16  eind  17  corresponding  to  case  14,  and  using  cutoff 
values  of  .90  and  .95  respectively,  and  again  in  Figures  18  and  19 
corresponding  to  case  5,  and  using  cutoff  values  of  .95  and  .975 
respectively.  This  implies  that  the  tracking  accuracy  increases  with 
a higher  cutoff  control  value,  for  the  environment  simulated  here,  and 
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that  even  correlation  coefficients  as  high  as  .95  do  not  guarantee 
accuracy.  The  direct  effect  of  raising  the  cutoff  control  value  is 
to  reject  the  data  in  more  rows  and  columns  specifically,  those  rows 
or  columns  that  did  not  correlate  sufficiently  high.  Hopefully,  this 
results  in  a more  careful  selection  of  "good"  data,  and  will  result 
in  better  tracking  accuracy.  Since  the  choice  of  a good  cutoff  value 
is  very  problem  dependent,  with  noisy  data  leading  to  lower  correlation 
coefficients,  this  value  v;ould  have  to  be  changed  on-line  for  the  best 


results . 
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VII . Conclusions  and  Recommendations 

Conclusions 

It  is  apparent  in  all  the  cases  studied  that  the  extended  Kalman 
filter  outperforms  the  correlation  tracker  under  the  same  set  of  input 
conditions.  As  shown  in  Tables  II  and  III,  the  performance  difference 
is  shown  most  dramatically  when  the  signal- to-noise  ratio  is  low,  and 
again  when  the  Gaussian  beam  dispersion  is  small  compared  to  the  pixel 
size.  This  result  was  expected,  however,  since  when  the  signal- to- 
noise  ratio  is  low,  there  is  little  definition  to  the  target  intensity 
pattern  and  the  tracking  algorithm  is  receiving  mostly  noisy  data. 

Similarly,  when  the  beam  dispersion  is  small,  all  of  the  useful  infor- 
mation is  concentrated  in  the  area  of  a few  pixels  and  is  easily  masked 
by  noise.  This  masking  is  especially  strong  at  low  signal-to-noise 
ratios . Since  the  correlation  algorithm  is  given  no  a priori  informa- 
tion about  either  the  motion  of  the  target  or  what  the  intensity 
pattern  looks  like,  it  does  a poor  job  of  picking  out  the  target. 

> 

Nevertheless,  this  type  of  tracker  could  be  effectively  used  in  an 
environment  where  the  signal  was  strong,  or  where  the  shape  of  the 
object  being  tracked  is  completely  unknown  or  constantly  changing. 

These  advantages  could,  in  fact,  lead  to  a design  that  is  even  more 
robust  than  the  extended  Kalman  filter. 

The  extended  Kalman  filter  performed  well  in  all  cases,  even  with 
a very  low  signal-to-noise  ratio.  Prior  knowledge  of  the  size,  shape, 
and  motion  of  the  target  was  used  to  track  the  target  effectively 
under  a variety  of  conditions.  The  filter  algorithm  also  appeared 
to  be  relatively  insensitive  to  small  changes  in  the  size  of  the 
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assumed  Gaussian  intensity  function.  T^iis  implies  that  the  filter 
design  may  be  rather  robust  and,  therefore,  could  be  used  in  an  en- 
vironment where  parameters  are  constantly  varying  or  are  not  well 
kncwn  a priori.  It  is  apparent,  however,  that  an  extremely  poor 
assumption,  such  as  assiming  an  entirely  false  shape,  could  signifi- 
cantly degrade  the  performance  of  this  algorithm.  This  implies  that 
this  type  of  filter  could  be  used  effectively  only  when  at  least  one 
aspect  of  the  target  were  known  with  some  degree  of  certainty. 

Recommendations 

Further  research  is  recommended  to  assess  the  robustness  of  the 
basic  design  of  the  filter ' presented  here.  These  tests  should  in- 
clude cases  where  there  are  errors  in  the  assumed  model  of  atmospheric 
turbulence  as  well  as  poor  assumptions  of  the  size,  shape,  and  motion 
of  the  target.  If  warranted,  this  analysis  could  lead  into  the  study 
or  design  of  an  adaptive  filter  model  that  could  estimate  parameters 
that  are  critical  to  the  filter's  performance  and  robustness,  such 
as  target  pattern  size  or  shape,  on-line.  It  could  even  be  used  to 
describe  the  motion  of  the  target  or  the  atmospheric  jitter  in  order 
to  adapt  to  changing  conditons  or  targets. 

Another  area  of  suggested  research  is  to  attempt  tracking  under 
a wide  variety  of  background  conditions.  In  this  way,  performance  of 
the  filter  could  be  examined  as  a function  of  background  clutter  as 
would  be  found  in  a "look  down"  mode. 

In  order  to  improve  the  computational  feasibility  of  the  filter 
for  on-line  applications,  research  should  be  performed  on  ways  to 
improve  the  speed  and  memory  requirements  of  this  type  of  algorithm. 
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One  way  to  improve  the  speed  would  be  to  incorporate  measurements 
using  only  a small  number  of  pixels . This  would  involve  detenning 
which  pixels  contain  the  most  information  so  o;ily  these  would  lie 
exploited.  One  v/ay  to  do  this  would  be  to  use  x(tT)  as  the  best 
guess  of  the  position  of  the  target  at  the  measurement  time  and  use 
only  those  pixels  that  lie  in  some  locus  centered  around  x.  Otlier 
methods , such  as  square  root  filtering  techniques  or  U-D  implementa- 
tion would  reduce  the  computer  wordlength  requirements  but  would  re- 
sult in  increased  memory  and  time  requirements , These  techniques 
would  probably  be  cf  little  benefit  here  since  the  extended  Kalman 
filter  only  has  four  states.  With  an  adaptive  filter,  however, 
these  techniques  may  be  helpful  since  several  more  parameters  would 
have  to  be  estimated. 

A last  possibility  would  be  to  explore  how  an  extended  Kalman 
filter  of  this  type  would  operate  when  there  is  more  than  one  target 
in  the  FLIR  field  of  view.  The  goal  here  would  be  to  see  if  the 
filter  could  adequately  distinguish  and  track  targets  in  a combat 


environment  with  a short  time  constraint. 
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Appendix  A 


Numerical  Precision  Problems 


In  the  early  stages  of  this  study,  a great  deal  of  time  was  spent 

searching  for  a reasonable  way  to  evaluate  the  truth  model  matrix, 

d 

defined  in  equation  (25)  , emd  repeated  here  as 


4+1 


2t. 


= / 


t. 
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£^(t.^l,T)Gg,(T)2,j,(T)  G^(T)^{t.^^,T)dT 


(A-ll 


Initially,  the  truth  model  state  equations  were  put  into  phase  vari- 
able form,  since  it  could  be  done  from  inspection  of  the  system  transfer 
functions . 

Because  the  truth  model  is  a time-invariant  system,  the  integral  in 
equation  (A-1)  can  be  solved  exactly,  and  will  be  a constant  for  all 
sampling  periods . This  means  that  it  need  be  calculated  only  once 
off-line.  Solving  equation  (A-1)  exactly,  however,  would  require  a 
significant  cimount  of  manual  calculation,  even  after  reducing  the 
problem  because  of  the  two-channel  symmetry.  In  order  to  save  both 
time  and  effort,  approximation  techniques  were  examined. 

The  simplest  approximation  would  be  to  let 


At 


(A-2) 


This  results  in  a good  approximation  when  the  sampling  period  is  small 
compared  to  the  response  time  of  the  system  (Ref  9 :Ch  6) . In  this  case, 
however,  using  the  phase  variable  form  of  the  state  equations,  where 
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t 


r 


0 Kab" 


0 0 0 


0 Kab^ 


this  approximation  results  in  a diagonal  matrix.  This  would  effectively 
eliminate  the  effects  that  would  normally  be  present  due  to  any  cross 
correlation  terms  in  the  fully  evaluated  the  approximation 

was,  therefore,  judged  inadequate. 

The  second  approximation , was  a trapezoidal  integration  of  equation 
(A-1)  yielding 


This  approximation  to  Qj  appeared  to  be  adequate  until  the  matrix  v;as 

d 

passed  through  the  Cholesky  square  root  algorithm,  as  described  in 
Chapter  II,  where  it  was  found  to  contain  a negative  eigenvalue.  Knowing 
that  the  matrix  should  always  be  positive  semi-definate  (Ref  9) , this 
was  also  found  to  be  an  inadequate  approximation. 

Oil  the  next  attempt,  the  equation 


2ir  2.T  Et 
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where  is  defined  in  (A-4) , was  integrated  using  an  Euler  integration 
with  ten  steps  per  time  period.  This  again  resulted  in  a matrix  with 
a negative  eigenvalue  when  passed  through  the  Cholesky  square  root 
algorithm.  Trying  to  gain  accuracy  by  using  twenty  steps  gave  almost 
identical  results  and  this  method  was  also  abandoned. 

At  this  point,  the  phase  variable  form  for  the  state  equations 
was  reexamined  to  check  the  validity,  and  then  transformed  into  the 
Jordan  canonical  form  to  see  if  better  results  could  be  obtained. 

The  approximations  in  equations  (A-5)  and  (A-6)  were  again  tried 
with  tlie  canonical  form,  shown  in  equation  (13) . Once  again  the  results 
were  similar. 

As  a last  effort,  the  exact  integration  of  equation  (A-1)  was 
attempted,  using  the  canonical  form  because  of  the  ease  in  calculating 
the  ir  (ti_^_i,ti)  matrix.  This  effort  was  successful,  and  the  results 
are  listed  in  the  program  listing  found  in  Appendix  D. 

There  is  no  apparent  reason  why  the  approximations  of  equation 
(A-3)  and  (A-4)  resulted  in  matrices  that  were  not  positive  semi- 
definate . The  60  bit  word  size  of  the  CDC-6600  computer  system  would 
normally  be  large  enough  to  compensate  for  any  round-off  errors  that 
would  occur.  This  probably  means  that  this  particular  system  is  ill- 
conditioned,  and  extremely  sensitive  to  minor  errors. 
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Appendix  B 


Trans fontiing  the  Atmospheric  Model  from  Phase 
Variable  to  Jordan  Canonical  Form 


The  third  order  atmoshperic  disturbance  model  as  shown  in  Figure 
5 has  a transfer  function  described  as 


G(s)  = 


Kab^ 

(s  + a)  (s  + b)  ^ 


Kab" 


s^  + (2b  + a)s^  + (b^  + 2ab)s  + ab^ 


(B-1) 


(B-2) 


Transferring  this  into  a vector  differential  equation  gives 


x(t) 


■ 0 

1 

0 

■ 0 ■ 

0 

0 

- (b^+2ab) 

1 

- (2b+a) 

x(t) 

+ 

0 

_Kab2_ 

w(t) 


(t)  ~ ° oj  }^{t) 


(B-3) 


(B-4) 


Which  is  tlie  phase  variable  form  of  the  differential  equation.  This 
equation  can  be  transfonaed  into  the  Jordan  canonical  form  through  the 
use  of  a similarity  transformation  (Ref  3:13)  . This  is  accomplished 
by  using  the  following  equations 


F*  = V"^F  V 


G*  = V 'g 


H*  = H V 


(B-5) 

(B-6) 

(B-7) 
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where  V is  a matrix  made  up  of  the  eigenvectors  cind  generalized 


eigenvector  corresponding  to  the  eigenvalues  of  the  third  order  system. 

In  this  case , the  eigenvalues  are  easily  definahle  from  the  trans- 
fer function  as  -a,  -b,  -b  . The  corresponding  eigenvectors 
and  the  generalized  eigenvector  £ can  then  be  determined  from  the 
following  equations : 


(F  + a£) 

= 0 

(B-8) 

(F  + b^)^ 

= 0 

(B-9) 

(F  + M)g 

-2 

(B-10) 

Solving  these  equations  for  the  eigenvalues  gives 

110 

-a  -b  1 (B-ll) 

b^  -2b 

and 

b^  2b  1 

(a^-2ab)  -2b  -1  1 (B-12) 

(a-b) ^ 

(a^b-ab^)  (a^-  b*)  (a-b) 


Preforirting  the  transformation  shewn  in  equations  (B-5  through  (B-7) 
gives 


G* 


Kab^ 

Ta-b)^ 

Kab^ 
(a-b) 2 

Kab^ 

(a-b) 


H*  = 


1 


o] 


(B-14) 


(B-15) 


The  transformed  vector  differential  equation  in  Jordan  canonical  form 


IS  now 


r Kab^ 

-a 

0 

0 

(a-b) 

x(t)  = 

0 

-b 

. 1 

x(t)  + 

Kab^ 
(a-b) 2 

0 

0 

-b 

Kab^ 

- 

(a-b) 

w(t) 


(B-16) 


y(t)  =•  [l  1 O]  x(t) 


(B-17) 
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Appendix  C 


Gain  Calculation  for  Atmospheric  Jitter  Model 


The  atmospheric  jitter  model  used  in  this  study  was  based  on  a 
report  written  by  the  Analytic  Sciences  Corporation  (TASC)  (Ref  1) . 
The  report  showed  that  the  jitter  could  be  adequately  described  as 
the  output  of  a third  order  filter  driven  by  zero  mean,  white, 
Gaussian  noise.  This  model,  as  presented  in  Figure  5,  is 


w 

Kab^ 

y 

(s+a) (s+b) ^ 

(C-1) 


where 

w = an  input  white  noise  process  with  zero  mean  and 
variance  equal  to  1; 

a,b  = break  frequencies  for  the  process  (a=14.14  rad/sec, 
b=659.5  rad/sec ; 

y = output  of  the  system; 

K = system  gain. 

Since  the  strength  of  the  white  noise  term  is  specified  as  1,  the  gain, 

K , must  be  adjusted  to  obtain  the  desired  RI-IS  jitter  characteristics 
on  the  output.  In  mathematical  form,  this  means  that  K is  adjusted  to 
obtain 


[y’(ti)] 


= 


(C-2) 


where  is  the  desired  RMS  jitter. 

The  Jordan  canonical  form  of  the  system  in  (C-1)  was  derived  in 


Appendix  B to  be 


1 


y(t)  = 

1 o]  x(t) 

(C-4) 

Using  this  form 

of  the  equation,  the  desired  RI^S  jitter  can 

now  be 

e:q)ressed  as 

- 

y(t)] 

(C-5) 

= 

(Xj^  (t)  + X2  (t) ) 

(C-6) 

= 

x|(t)  + 2Xj^(t)x2(t)  + 

(C-7) 

= 

^ 1 

x|(t)J  + 2E  ^X]^(t)x2(t)J  + E [^X2(t)J 

(C-8) 

The  variance 

of 

the  desired  jitter  characteristics  is  equal 

to  the  sum 

of  the  variances  of  the  first  two  states  in  equation  (C-3)  plus  twice 
the  covariance  between  the  states,  all  at  time  t . 

The  continuous  time  model  for  the  covariance  matrix  can  be  written 
as 

P^t)  = FP(t)  + P(t)F^  + (C-9) 

for  this  case  where 

P(t)  = the  covariance  matrix 
F = the  system  plant  matrix 

-r 

t £ = the  system  input  matrix 

Q = the  variance  of  the  white  noise  process  driving 
the  system 
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when  the  system  has  reached  steady  state,  the  P(t)  matrix  will  equal 


zero,  and  equation  (C-9)  can 

be  written  as 

O = FP(t)  + 

P(t)F^  + GQg'^ 

(C-10) 

Substituting  the  appropriate 

F , G , and  Q 

into  equation  (C- 

-10) 

yields  the  equation 

-2aP-j^]^+c 

(s+t))  +E 

o = 

-P2l(a+b)+Pi3-c 

-2bP22+P32+P23‘'’'^ 

-2bP23+P33-E 

(c-11) 

-Pjj^  (a+b)+E 

-2bP32+P33-E 

-2bP33+D 

where 


c = K^a^b"*/ (a-b)  ** 


E = K^a^bVCa-b)^ 

D = K^a^bVCa-b)^ 

Making  use  of  the  fact  that  the  covariance  matrix  is  always  symmetric, 
equation  (C-11)  can  be  solved  to  give  the  steady  state  values  of  the 
covariance  matrix.  Solving  for  the  values  of  interest  yields 

K^ab-* 


E [Xj^ct)]  = = 


2 (a-b) 


E [xi(t)x2(t)]  = P^2  " ' 


(a+b)  (a-b) 


(a+b)  (a-b) 


(C-12) 

(C-13) 


E [^(t)]  = P 


22 


„2  2.  „2  2,  2 ^2  2,  3 

Kab  Kab  Kab 

'3  + 


4(a-b)^  2(a-b)^  2 (a-b)  '* 


(C-14) 


These  values  can  now  be  substituted  into  equation  (C-8)  along  with  the 
true  values  of  a eind  b to  give 
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K .382109544 


» 

(C-15) 


This  equation  will  result  in  the  correct  value  for  gain  given  a desired 
RMS  jitter  in  units  of  pixels. 
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Program  Listing 

This  appendix  contains  the  FORTRAN  source  code  for  the  computer 
program  used  in  this  study . The  program  was  written  for  implementation 
on  a CDC-6600  computer  system  which  has  a 60  bit  wordlength. 
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Appendix  E 


Plotted  Outputs  of  Cases  Studied 

This  appendix  contains  the  plotted  outputs  for  the  sixteen  cases 
studied.  There  are  three  plots  for  each  case; 

- mean  error  ±1  sigma  for  the  correlator 

- mean  error  ±1  sigma  for  the  filter  estimate  of  target 
dynamics 

- mean  error  ±1  sigma  for  the  filter  estimate  of  atmospheric 
jitter 

Because  of  the  similarity  in  results  between  the  horizontal  (x)  and 
vertical  (y)  channels , only  the  mean  error  plots  for  the  horizontal 
channel  are  presented. 
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Figure  20a. 
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Figure  20c.  X CHANNEL  ATMOSPHERICS  ERROR  (S/N=20] 


I 


MEPN  FROCK I NG  ERROR 
FOR  CORRELATOR  MODEL 


09*0 


00*0 

S13XId 


09*0- 


00*  I 


00*T- 


nCPN  ERROR  4-1  SlGnP 
PTMOSPMERIC  OYNPnlCS 
X CMPNNEL 
FLIR  fOV 


09  ^ 

08*0  00 '0 

08*6- 

^<J 

09’T- 

S13XId 

Figure  22c.  X CHANNEL  ATMOSPHERICS  ERROR  (5/N=20) 
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Figure  25c.  X CHANNEL  ATMOSPHERICS  ERROR  (S/N=10) 
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Figure  26a.  X CHRNNEL  CORRELATOR  TRACKING  ERROR  ( S/N  = 10) 
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Figure  26b.  X CHANNEL  OYNflMICS  ERROR  (S/N=lO) 
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Figure  27a.  X CHANNEL  CORRELATOR  TRACKING  ERROR  ( 5/N  = 10  J 
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Figure  27b.  X CHRNNEL  DYNAMICS  ERROR  (S/N=10) 
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Figure  31b.  x CHANNEL  DYNAMICS  ERROR  (5/N=  1 ) 
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